function [PartProb,Npdf,obsNv,maxP_TS,EProf,EN,Prof] = Full_Profit(muT,sigmaT,sigmaF0,sigmaF1,rho,mu_F_P,mu_F_TS,A,R,TG,rp,snPass,iobs)

N       = 0:1:50;
Nbar    = 35;
a       = 1.38;
b       = 0.82342;

if rp.parl==0
    [s1,s2,s3,s4] = randdiv2(snPass.s,iobs,rp.k,rp.m);
    
    N3D = snPass.N3D;
    Nv  = snPass.Nv;
    %Ngeq3D = snPass.Ngeq3D;
    N03D = snPass.N03D;
    %Ngeqv = snPass.Ngeqv;
    N0v   = snPass.N0v;
elseif rp.parl==1
    rng(rp.theseed+iobs);
    s = rand(3*rp.k + rp.m,1);

    s1 = s(1:rp.k,1);              % Size k for T at t=0  
    s2 = s(rp.k + 1:rp.k + rp.m,1)';         % Size m for F1 at t=0
    s2 = s2(ones(1,rp.k),:);       % Repeated k times for a kxm vector
    s3 = s(rp.k + rp.m + 1:2*rp.k + rp.m,1);      % Size k for F0 at t=0  
    s4 = s(2*rp.k + rp.m + 1:3*rp.k + rp.m,1);    % Size k for F1 at t=1
    
    Ntr     = permute(N,[3,1,2]);
    N3D     = Ntr(ones(rp.k,1),ones(rp.m,1),:);
    Nv      = N(ones(rp.k,1),:);
    %Ngeq3D  = (N3D>=Nbar);
    N03D    = (N3D>0);
    %Ngeqv   = (Nv>=Nbar);
    N0v     = (Nv>0);
end

% For stochastic calcs ***
Nst     = zeros(0,51);
for ii=N
    Nst(ii+1) = 1-bbcdf(ii,Nbar-1,a,b);
end
Nsttr   = permute(Nst,[3,1,2]);
Ngeq3D  = Nsttr(ones(rp.k,1),ones(rp.m,1),:);
Ngeqv   = Nst(ones(rp.k,1),:);
% ************************


etaT = norminv(s1,0,1);
eta1 = norminv(s2,0,1);     % Size k x m ~ unknown farmer at time 0
eta2 = norminv(s3,0,1);     % Length k ~ unknown to the econometrician
eta3 = norminv(s4,0,1);     % Length k ~ unknown to the econometrician

s21  = rho.*sigmaF0;
s22  = sqrt(sigmaF0^2.*(1-rho^2));

T     = exp(muT+etaT.*sigmaT);
F0_P  = etaT.*s21 + eta3.*s22 + mu_F_P;  % Used for participation dec
F0_TS = etaT.*s21 + eta3.*s22 + mu_F_TS; % Used for tree number dec
F1ind = eta1.*sigmaF1;       % Used for participation dec
F1    = eta2.*sigmaF1;       % Used for tree number dec

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%  Produce 3D comformable matrices so that we can  ;
%  can conduct the optimization in a single step   ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;

% T3D_P is a (k x m x 51) matrix that has 5 draws of T repeated
% (m x 51) times;
T3D_P  = T(:,ones(rp.m,1),ones(51,1));

% F03D_P is a (k x m x 51) matrix that has k draws of F0 repeated
% (m X 51) times;
F03D_P = F0_P(:,ones(rp.m,1),ones(51,1));

% F13D_P is a (k x m x 51) matrix that has m draws of F1 repeated
% (m x 51) times;
F13D_P = F1ind(:,:,ones(51,1));

% T_TS, F1_TS and F0_TS are T, F0 and F1 vectors of k draws
% repeated 51 times;
T_TS   = T(:,ones(1,51));
F0_TS  = F0_TS(:,ones(1,51));
F1_TS  = F1(:,ones(1,51));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
% PARTICIPATION DECISION                           ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;

if TG==1
    Prof = N3D - (1./(2*T3D_P)).*(N3D.^2) - N03D.*(F03D_P+F13D_P);
else
    Prof = N3D - (1./(2*T3D_P)).*(N3D.^2) - N03D.*(F03D_P+F13D_P) ...
        + Ngeq3D.*R;
end

%--------------------------------------------------;
% Generate maximum profit and optimal choice of    ;
% trees                                            ;
%--------------------------------------------------;

[maxP_P,argmaxN] = max(Prof,[],3);
obsN = argmaxN - 1;

%-------------------------------------------------;
% Participation, Expected Profit and Survival      ;
% outcomes.                                        ;
%--------------------------------------------------;

EProf = mean(maxP_P,2);
EN    = mean(obsN,2);

Part  = (rp.df.*EProf>=12-A);   % Participation rates;

PartProb = mean(Part,1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
% TREE SURVIVAL DECISION                           ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;

%---------------------------------------------------;
% pdf for number of trees                           ;
%---------------------------------------------------;

Profv = Nv - (1./(T_TS.*2)).*(Nv.^2) - N0v.*(F0_TS+F1_TS) ...
    + Ngeqv.*R;

Profv = Profv(Part,:);

if isempty(Profv)==1
    maxP_TS  = 0;
    obsNv = 0;
    Npdf     = [1 zeros(1,50)];
else
    if rp.lam>0
        maxP_TS = max(Profv,[],2);
        Profvt  = bsxfun(@minus,Profv,maxP_TS); % Normalization
        expProf = exp(Profvt./rp.lam);    % Constructing continuous functions for the probability of
        denom   = sum(expProf,2);               % choosing each value of N;
        denom   = denom(:,ones(1,51));
        logitProf = expProf./denom;             % Computes the probability of choosing each value of N;
        Npdf    = mean(logitProf,1);         % We average across draws.
    else
        [maxP_TS,argmaxNv] = max(Profv,[],2);
        obsNv   = binornd(argmaxNv - 1,betarnd(a,b));
        Npdf    = mean(bsxfun(@eq,obsNv,N),1);
    end
end

